Here are the packages I will need:

library(ggplot2)
library(dplyr)
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(tidyr)
library(tidyverse)
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ forcats   1.0.0     ✔ readr     2.1.4
## ✔ lubridate 1.9.3     ✔ stringr   1.5.1
## ✔ purrr     1.0.2     ✔ tibble    3.2.1
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(skimr)
library(cowplot)
## 
## Attaching package: 'cowplot'
## 
## The following object is masked from 'package:lubridate':
## 
##     stamp

#Loading and reading data

#Challenge 1

One-Factor ANOVA and Inference Step 1 Make boxplots of log(Mass) in relation to Trophic.Level and Migration behavior type. For each plot, drop from the visualization all species records where the categorical variable of interest is missing from the dataset. Also, you will want to convert the variable Migration (which is scored as a number: “1”, “2”, or “3”) from class numeric to either being classified as a factor or as a character (string) variable.

#Make boxplots of log(Mass) in relation to Trophic.Level and Migration behavior type.
ggplot(data=d |> drop_na(Trophic.Level), aes(x = Trophic.Level, y=log(Mass)))+ geom_boxplot() + geom_jitter( alpha = 0.05)

ggplot(data=d|> drop_na(Migration), aes(x=as.factor(Migration), y=log(Mass)))+geom_boxplot()

#an alternate way below
d <- d |>
    mutate(logMass = log(Mass), logRS = log(Range.Size), logBeak = log(Beak.Length_Culmen),
        logTarsus = log(Tarsus.Length), Migration = as.factor(Migration))

p1 <- ggplot(data = d |>
    drop_na(Trophic.Level), aes(x = Trophic.Level, y = log(Mass))) + geom_boxplot() +
    geom_jitter(alpha = 0.05)
plot(p1)

p2 <- ggplot(data = d |>
    drop_na(Migration), aes(x = Migration, y = log(Mass))) + geom_boxplot() + geom_jitter(alpha = 0.05)

#plotting the box plots next to eachother
plot_grid(p1, p2, nrow = 1)

Step 2 Run linear models using the lm() function to look at the relationship between log(Mass) and Trophic.Level and between log(Mass) and Migration.

Is log(Mass) associated with either Trophic.Level or Migration category? That is, in the global test of significance, is the F statistic large enough to reject the null hypothesis of an F value of zero? Yes, the global test is indeed significant for both Trophi.Level and Migration.

Given the regression coefficients returned for your Migration model, which Migration categor(ies) are different than the reference level? Levels 2 and 3 differ from the reference level. What level is the reference level? The reference level is level 1. Relevel and assess differences among the remaining pair of Migration categories. Levels 2 and 3 differe again after releveling.

m1 <- lm(log (Mass) ~ Trophic.Level, data=d)

m2 <- lm(log(Mass) ~ as.factor(Migration), data = d)

m1
## 
## Call:
## lm(formula = log(Mass) ~ Trophic.Level, data = d)
## 
## Coefficients:
##            (Intercept)  Trophic.LevelHerbivore   Trophic.LevelOmnivore  
##                3.80834                 0.25639                 0.01422  
## Trophic.LevelScavenger  
##                4.63189
m2
## 
## Call:
## lm(formula = log(Mass) ~ as.factor(Migration), data = d)
## 
## Coefficients:
##           (Intercept)  as.factor(Migration)2  as.factor(Migration)3  
##                3.7746                 0.7597                 0.3765
summary(m1)
## 
## Call:
## lm(formula = log(Mass) ~ Trophic.Level, data = d)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.4229 -1.1551 -0.3028  0.8982  7.5526 
## 
## Coefficients:
##                        Estimate Std. Error t value Pr(>|t|)    
## (Intercept)             3.80834    0.01967 193.632  < 2e-16 ***
## Trophic.LevelHerbivore  0.25639    0.03406   7.528 5.54e-14 ***
## Trophic.LevelOmnivore   0.01422    0.04116   0.345     0.73    
## Trophic.LevelScavenger  4.63189    0.34447  13.446  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.538 on 11000 degrees of freedom
##   (5 observations deleted due to missingness)
## Multiple R-squared:  0.02094,    Adjusted R-squared:  0.02067 
## F-statistic: 78.42 on 3 and 11000 DF,  p-value: < 2.2e-16
summary(m2)
## 
## Call:
## lm(formula = log(Mass) ~ as.factor(Migration), data = d)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.8924 -1.1769 -0.3088  0.9152  7.8427 
## 
## Coefficients:
##                       Estimate Std. Error t value Pr(>|t|)    
## (Intercept)            3.77457    0.01636 230.710  < 2e-16 ***
## as.factor(Migration)2  0.75971    0.04731  16.059  < 2e-16 ***
## as.factor(Migration)3  0.37647    0.05155   7.303 3.02e-13 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.535 on 10983 degrees of freedom
##   (23 observations deleted due to missingness)
## Multiple R-squared:  0.02563,    Adjusted R-squared:  0.02546 
## F-statistic: 144.5 on 2 and 10983 DF,  p-value: < 2.2e-16
#Relevel
#MMigration <- as.factor(d$Migration)
#d <- d|> mutate(Migration = relevel (MMigration, ref = "3"))
d <- d|> mutate(Migration = relevel (Migration, ref = "3"))
m2 <- lm(log(Mass) ~ Migration, data = d)

summary(m2)
## 
## Call:
## lm(formula = log(Mass) ~ Migration, data = d)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.8924 -1.1769 -0.3088  0.9152  7.8427 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  4.15104    0.04889  84.909  < 2e-16 ***
## Migration1  -0.37647    0.05155  -7.303 3.02e-13 ***
## Migration2   0.38324    0.06603   5.804 6.67e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.535 on 10983 degrees of freedom
##   (23 observations deleted due to missingness)
## Multiple R-squared:  0.02563,    Adjusted R-squared:  0.02546 
## F-statistic: 144.5 on 2 and 10983 DF,  p-value: < 2.2e-16

Step 3 Conduct a post-hoc Tukey Honest Significant Differences test to also evaluate which Migration categories differ “significantly” from one another (see Module 20).

m1aov <- aov(log(Mass) ~Trophic.Level, data=d)
(pairwise.t.test(log(d$Mass), d$Trophic.Level, p.adj = "bonferroni"))
## 
##  Pairwise comparisons using t tests with pooled SD 
## 
## data:  log(d$Mass) and d$Trophic.Level 
## 
##           Carnivore Herbivore Omnivore
## Herbivore 3.3e-13   -         -       
## Omnivore  1         6.7e-07   -       
## Scavenger < 2e-16   < 2e-16   < 2e-16 
## 
## P value adjustment method: bonferroni
(posthoc <- TukeyHSD (m1aov, which= "Trophic.Level", conf.level = 0.95))
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = log(Mass) ~ Trophic.Level, data = d)
## 
## $Trophic.Level
##                            diff         lwr        upr     p adj
## Herbivore-Carnivore  0.25638774  0.16888205  0.3438934 0.0000000
## Omnivore-Carnivore   0.01422185 -0.09154436  0.1199881 0.9858317
## Scavenger-Carnivore  4.63188503  3.74679736  5.5169727 0.0000000
## Omnivore-Herbivore  -0.24216589 -0.35936711 -0.1249647 0.0000007
## Scavenger-Herbivore  4.37549729  3.48897045  5.2620241 0.0000000
## Scavenger-Omnivore   4.61766318  3.72914808  5.5061783 0.0000000
#(posthoc <- TukeyHSD (m1aov, which= "Migration", conf.level = 0.95))
plot(posthoc)

m2 <- aov(log(Mass) ~ Migration, data = d)
(posthoc <- TukeyHSD(m2, which = "Migration", conf.level = 0.95))
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = log(Mass) ~ Migration, data = d)
## 
## $Migration
##           diff        lwr        upr p adj
## 1-3 -0.3764693 -0.4973105 -0.2556282     0
## 2-3  0.3832374  0.2284536  0.5380211     0
## 2-1  0.7597067  0.6488157  0.8705977     0
plot(posthoc)

Step 4

Use a permutation approach to inference to generate a null distribution of F statistic values for the model of log(Mass) in relation to Trophic.Level and calculate a p value for your original F statistic. You can do this either by programming your own permutation test (e.g., by shuffling values for the predictor or response variable among observations and calculating an F statistic for each replicate) or by using the {infer} workflow and setting calculate(stat=“F”).

library(broom)
original.F <- aov(log(Mass) ~Trophic.Level, data=d) |>
  tidy() |>
  filter(term== "Trophic.Level")|>
  pull(statistic)

original.F
## [1] 78.42283
#by using infer
library(infer)

d <-d|> mutate(logMass = log(Mass))
permuted.F <-d|>
  specify(logMass ~Trophic.Level)|>
  hypothesize(null = "independence") |>
  generate(reps = 1000, type = "permute") |>
  calculate(stat ="F")
## Warning: Removed 5 rows containing missing values.
visualize(permuted.F) + shade_p_value(obs_stat = original.F, direction = "greater")
## Warning in (function (mapping = NULL, data = NULL, stat = "identity", position = "identity", : All aesthetics have length 1, but the data has 1000 rows.
## ℹ Did you mean to use `annotate()`?
## Warning in min(diff(unique_loc)): no non-missing arguments to min; returning Inf

p_value<- permuted.F|> get_p_value(obs_stat = original.F, direction = "greater")
## Warning: Please be cautious in reporting a p-value of 0. This result is an approximation
## based on the number of `reps` chosen in the `generate()` step.
## ℹ See `get_p_value()` (`?infer::get_p_value()`) for more information.
#original.F$p.value

#permuted.F <- as_tibble(permuted.F) |>
 #   rename(stat = "result")

p.value <- permuted.F |>
    get_p_value(obs_stat = original.F, direction = "greater")
## Warning: Please be cautious in reporting a p-value of 0. This result is an approximation
## based on the number of `reps` chosen in the `generate()` step.
## ℹ See `get_p_value()` (`?infer::get_p_value()`) for more information.
p.value
## # A tibble: 1 × 1
##   p_value
##     <dbl>
## 1       0

Challenge 2

Data Wrangling, One- and Two-Factor ANOVA

Step 1 Create the following two new variables and add them to AVONET dataset: Relative beak length, which you should calculate as the residual of log(Beak.Length_Culmen) on log(Mass). Relative tarsus length, which you should calculate as the residual of log(Tarsus.Length) on log(Mass).

d$relBL <-resid(lm(formula=log(Beak.Length_Culmen)~log(Mass), data=d))

d$relTL <-resid(lm(log(Tarsus.Length)~log(Mass),data=d))

relBeak <- d$relBL

relTarsus <- d$relTL

Step 2 Make a boxplot or violin plot of your new relative tarsus length variable in relation to Primary.Lifestyle and of your new relative beak length variable in relation to Trophic.Niche

ggplot(d|>filter(!is.na(Primary.Lifestyle)),aes(x=Primary.Lifestyle, y=relTL))+
  geom_boxplot()+
  theme(axis.text.x=element_text(angle=90))

ggplot(data=d |> filter (!is.na(Trophic.Niche)), aes(x=Trophic.Niche, y = relBL))+
  geom_boxplot()+
  theme(axis.text.x=element_text(angle=90))

Step 3 Run an ANOVA analyses to look at the association between geographic range size and the variable Migration. You should first drop those observations for which Migration is not scored and also look at the distribution of the variable Range.Size to decide whether and how it might need to be transformed. Based on the global model, is range size associated with form of migration? How much of the variance in your measure of range size is associated with Migration behavior style?

Given the regression coefficients returned in output of the model, which Migration categor(ies) are different than the reference level? What level is the reference level? Relevel and assess differences among the remaining pair of Migration categories. Also conduct a post-hoc Tukey Honest Significant Differences test to also evaluate which Migration categories differ “significantly” from one another (see Module 20).

ggplot (d, aes(x=as.factor(Migration),y=Range.Size))+
  geom_boxplot()+geom_jitter()
## Warning: Removed 57 rows containing non-finite outside the scale range
## (`stat_boxplot()`).
## Warning: Removed 57 rows containing missing values or values outside the scale range
## (`geom_point()`).

m<-aov(Range.Size~as.factor(Migration), data=d)

summary(m)
##                         Df    Sum Sq   Mean Sq F value Pr(>F)    
## as.factor(Migration)     2 5.331e+16 2.666e+16   499.1 <2e-16 ***
## Residuals            10934 5.840e+17 5.341e+13                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 72 observations deleted due to missingness
TukeyHSD(m)
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = Range.Size ~ as.factor(Migration), data = d)
## 
## $`as.factor(Migration)`
##           diff      lwr        upr     p adj
## 1-3 -5948065.8 -6525065 -5371066.7 0.0000000
## 2-3  -780559.2 -1519665   -41453.6 0.0355345
## 2-1  5167506.6  4638082  5696930.7 0.0000000
ggplot(d|>filter(!is.na(Migration)), aes(x=as.factor(Migration), y=log(Range.Size)))+
  geom_boxplot() +geom_jitter(aes(alpha=0.005))
## Warning: Removed 49 rows containing non-finite outside the scale range
## (`stat_boxplot()`).
## Warning: Removed 49 rows containing missing values or values outside the scale range
## (`geom_point()`).

m1reg <- lm(log(Range.Size) ~ as.factor (Migration),data=d)

summary(m1reg)|> tidy()
## # A tibble: 3 × 5
##   term                  estimate std.error statistic   p.value
##   <chr>                    <dbl>     <dbl>     <dbl>     <dbl>
## 1 (Intercept)             14.6      0.0890    164.   0        
## 2 as.factor(Migration)1   -2.52     0.0938    -26.8  1.14e-153
## 3 as.factor(Migration)2   -0.732    0.120      -6.10 1.13e-  9
#m1<-aov(log(Range.Size)~as.factor(Migration), data=d)

#summary(m1)
#migration
#Sedentary = 1
#partially migratory=2
#migratopry=3

Step 4 Winnow your original data to just consider birds from the Infraorder “Passeriformes” (song birds).

Run separate one-factor ANOVA analyses to look at the association between [1] relative beak length and Primary.Lifestyle and between [2] relative beak length and Trophic.Level. In doing so…

Make boxplots of response variable by each predictor and by the combination of predictors.

Run linear models for each predictor separately and interpret the model output. Step 5 Run a two-factor model to look at the association between relative beak length and both Primary.Lifestyle and Trophic.Level among the passeriforms. Based on the model output, what would you conclude about how relative beak length is related to these two variables? Step 6 Finally, run an additional two-way model with the same dataset and predictors, but adding the possibility of an interaction term. To do this, you should modify your model formula using the colon operator (:) to specify the interaction, e.g., relative beak length ~ Primary.Lifestyle + Trophic.Level + Primary.Lifestyle:Trophic.Level. Based on the model output, what would you now conclude about how relative beak length is related to these two variables? Step 7 Use the interaction.plot() function to visualize the interaction between Primary.Lifestyle and Trophic.Level (see Module 20).

p <-d|> filter(Order1 == "Passeriformes")

ggplot(d|>filter(!is.na(Primary.Lifestyle)), aes(x=Primary.Lifestyle, y= relBL))+
  geom_boxplot()

ggplot(d|>filter(!is.na(Trophic.Level)), aes(x=Trophic.Level, y= relBL))+
  geom_boxplot()

m0 <- aov(relBL~1, data=p)

m1<- aov(relBL ~ Primary.Lifestyle, data = p)
m2<- aov(relBL ~ Trophic.Level, data = p)
m3<- aov(relBL ~ Trophic.Level +Primary.Lifestyle + Trophic.Level:Primary.Lifestyle, data = p)
summary(m)
##                         Df    Sum Sq   Mean Sq F value Pr(>F)    
## as.factor(Migration)     2 5.331e+16 2.666e+16   499.1 <2e-16 ***
## Residuals            10934 5.840e+17 5.341e+13                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 72 observations deleted due to missingness
####
ggplot(data=p, aes(x=Primary.Lifestyle, y= relBL))+
  geom_boxplot()+
  geom_jitter(alpha =0.05)

ggplot(data=p, aes(x=Trophic.Level, y= relBL))+
  geom_boxplot()+
  geom_jitter(alpha =0.05)

#or we could facet wrap

ggplot(data=p, aes(x=Primary.Lifestyle, y= relBL))+
  geom_boxplot()+
  facet_wrap(vars(Trophic.Level))+
  geom_jitter(alpha =0.05)

ggplot(data=p, aes(x=Trophic.Level, y= relBL))+
  geom_boxplot()+
  facet_wrap(vars(Primary.Lifestyle))+
  geom_jitter(alpha =0.05)

interaction.plot(
  x.factor=p$Trophic.Level,
  xlab="Trophic Level",
  trace.factor = p$Primary.Lifestyle,
  trace.label="Primary.Lifestyle",
  response=p$relBL,
  ylab="Mean Relative Beak Length"
)

interaction.plot(
  x.factor=p$Primary.Lifestyle,
  xlab="PrimaryLifestyle",
  trace.factor = p$Trophic.Level,
  trace.label="Trophic Level",
  response=p$relBL,
  ylab="Mean Relative Beak Length"
)

m3 <-aov(relBL~Primary.Lifestyle+Trophic.Level, data=p)
m4<-aov(relBL~Primary.Lifestyle+Trophic.Level+ Primary.Lifestyle:Trophic.Level, data=p)

anova(m0,m1, test="F")
## Analysis of Variance Table
## 
## Model 1: relBL ~ 1
## Model 2: relBL ~ Primary.Lifestyle
##   Res.Df    RSS Df Sum of Sq      F    Pr(>F)    
## 1   6613 326.12                                  
## 2   6610 307.92  3      18.2 130.23 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(m0,m2, test="F")
## Analysis of Variance Table
## 
## Model 1: relBL ~ 1
## Model 2: relBL ~ Trophic.Level
##   Res.Df    RSS Df Sum of Sq      F    Pr(>F)    
## 1   6613 326.12                                  
## 2   6611 309.81  2    16.309 174.01 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(m1,m3, test="F")
## Analysis of Variance Table
## 
## Model 1: relBL ~ Primary.Lifestyle
## Model 2: relBL ~ Primary.Lifestyle + Trophic.Level
##   Res.Df    RSS Df Sum of Sq      F    Pr(>F)    
## 1   6610 307.92                                  
## 2   6608 290.24  2    17.677 201.23 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(m2,m3, test="F")
## Analysis of Variance Table
## 
## Model 1: relBL ~ Trophic.Level
## Model 2: relBL ~ Primary.Lifestyle + Trophic.Level
##   Res.Df    RSS Df Sum of Sq     F    Pr(>F)    
## 1   6611 309.81                                 
## 2   6608 290.24  3    19.568 148.5 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Below are incorporated class notes

fIn class #```{r}

d<- d |> mutate(logMass = log(Mass), logRS = log(Range.Size), logBeak = log(Beak.Length_Culmen), logTarsus= log(Tarsus.Length), Migration = as.factor(Migration))

m1 <-lm(data=d, logBeak ~logRS * Migration)

m2<-lm(data=d, logBeak ~logRS + Migration)

m3 <-lm(data=d, logBeak ~logRS)

m4 <- lm(data=d, logBeak ~Migration) m5 <- lm(data=d, logBeak ~1)

#anova(m2, m1, test = “F”)

#models didnt converge because one of them have data for

#we are comparting rs and migration vs rs anova(m3, m2, test = “F”)

#adding migration does the power.


With same AVONET dataset we will Explore forward and backward selection

#```{r}

d_new <- d %>%
  drop_na(logBeak, logRS, Migration, Trophic.Level, logTarsus, Primary.Lifestyle)
  
lm(data=d_new, logBeak ~ logRS + Migration +Trophic.Level + relTarsus + Primary.Lifestyle)

m_null <- lm(data=d_new, logBeak ~1)

add1(m_null, scope = .~. +logRS + Migration + logTarsus +Trophic.Level + Primary.Lifestyle, test = "F")

#Add primary lifestyle as the new model, it's the largest one

m1 <- update(m_null, formula = .~. + logTarsus)

summary(m1)
#low p value
add1(m1, scope = .~. +logRS + Migration +Trophic.Level + logTarsus + Primary.Lifestyle, test = "F")

m2<-update(m1, formula = .~. + Primary.Lifestyle)
summary(m2)
#overall still significant model multiple r2 inscreading

add1(m2, scope = .~. +logRS + Migration +Trophic.Level + logTarsus + Primary.Lifestyle, test = "F")

#trophic highest f value, 
m3<-update(m2, formula = .~. + Trophoic.Level)
m3
summary(m3)
#still significant cuz its getting larger


#intercept only model, which predictor adds the most explanitory value. until adding new predictors does not add more explanitory value
add1(m3, scope = .~. +logRS + Migration +Trophic.Level +logTarsus + Primary.Lifestyle, test = "F")

#trophic highest f value, 
m4<-update(m3, formula = .~. + logRS)
summary(m4)

Additional Steps? In the exercise above, we really did not do any checking with this dataset to see if the data meet the primary assumptions for standard linear regression and ANOVA, which are that variables/residuals within each grouping level are roughly normally distributed and have roughly equal variances. Sample sizes within each grouping level should also be roughly equal. As noted in Module 20, a general rule of thumb for “equal” variances is to compare the largest and smallest within-grouping level standard deviations and, if this value is less than 2, then it is often reasonable to presume the assumption may not be violated.

Use this approach to see whether variances in across groups in your various models (e.g., for relative beak length ~ trophic level) are roughly equal. Additionally, do a visual check of whether observations and model residuals within groups look to be normally distributed.

In class Backwards selection: ####```{r}

m_full <-lm(data=d_new, logBeak ~1 + logRS +logTarsus + Migration +Trophic.Level + Primary.Lifestyle)

drop1(m_full, test= “F”)

#since migration is the lowest F

m2 <-update(m_full, .~. -Migration) summary(m2)

drop1(m2, test= “F”)

m3 <-update(m2, .~. - logRS)

summary(m3)

drop1(m3, test= “F”)

m4 <- update(m3, .~. )

library(MASS)

highest AIC = best data

April 3

###```{r}
d_new <- d %>%
  drop_na(logBeak, logRS, Migration, Trophic.Level, logTarsus, Primary.Lifestyle)

m_full <- lm(data= d_new, relBeak ~ logRS + relTarsus + Migration + Trophic.Level + Primary.Lifestyle)


s <- stepAIC(m_full, scope=.~., direction = "both", trace = TRUE)

s<- stepAIC(m_null, scope = .~.)

####```{r}

library(MuMIn) m_full <- lm(data= d_new, relBeak ~ logRS + Migration + Trophic.Level + relTarsus + Primary.Lifestyle, na.action = na.fail) mods <- dredge(m_full) class(mods) mods.res <- get.models(mods, subset=delta <=4) #returns top models where delta.aicc <=4

mods.avg<-summary(model.avg(mods, subset=delta<= 4, fit = TRUE)) #returns top

mods.avg<-summary(model.avg(mods, subset = cumsum)) #mods.res$31

confint(mods.avg) plot(mods.avg, full = TRUE) plot(mods.avg, full = FALSE)

mods.avg <- summary (model.avg(mods, subset=delta <=4, fit = TRUE)) mods





####```{r}
f <- "https://raw.githubusercontent.com/difiore/ada-2024-datasets/main/Mammal_lifehistories_v2.txt"
d<- read_tsv(f, col_names= TRUE)

d[d==-999.00] <- NA
d<- dplyr::select(d, -`refs`)
d<- dplyr::select(d, -`litter size`)

colnames(d)
cols <- c("mass(g)", "gestation(mo)", "newborn(g)","weaning(mo)", "wean mass(g)","AFR(mo)","max. life(mo)","litters/year")

d <- d %>% mutate(across(all_of(cols), log))
 #replace 999
# drop refs and litter size variables
#log transform all other variables
#regress gestation, weaning, age at firs reproduction, and max lifespan on mass and add residuals to the dataframe (hint na.action=na.exclude)
#plot residuals of max lifespan in relation to order - which orders have the highest relative newborn mass

#which order have the highest residual lifespan